Genetic diversity and recombination of bovine enterovirus strains in China

ABSTRACT Bovine enterovirus (BEV) consisting of enterovirus species E (EV-E) and F (EV-F) is the causative agent associated with respiratory and gastrointestinal diseases in cattle. Here, we reported the characterization, genetic diversity, and recombination of novel BEV strains isolated from the major cattle-raising regions in China during 2012–2018. Twenty-seven BEV strains were successfully isolated and characterized. Molecular characterization demonstrated that the majority of these novel BEV strains (24/27) were EV-E, while only few strains (3/27) were EV-F. Sequence analysis revealed the diversity of the circulating BEV strains such as species and subtypes where different species or subtype coinfections were detected in the same regions and even in the same cattle herds. For the EV-E, two novel subtypes, designated as EV-E6 and EV-E7, were revealed in addition to the currently reported EV-E1–EV-E5. Comparative genomic analysis revealed the intraspecies and interspecies genetic exchanges among BEV isolates. The representative strain HeN-B62 was probably from AN12 (EV-F7) and PS-87-Belfast (EV-F3) strains. The interspecies recombination between EV-E and EV-F was also discovered, where the EV-F7-AN12 might be from EV-E5 and EV-F1, and EV-E5-MexKSU/5 may be recombined from EV-F7 and EV-E1. The aforementioned results revealed the genetic diversity and recombination of novel BEV strains and unveiled the different BEV species or subtype infections in the same cattle herd, which will broaden the understanding of enterovirus genetic diversity, recombination, pathogenesis, and prevention of disease outbreaks. IMPORTANCE Bovine enterovirus (BEV) infection is an emerging disease in China that is characterized by digestive, respiratory, and reproductive disorders. In this study, we first reported two novel EV-E subtypes detected in cattle herds in China, unveiled the coinfection of two enterovirus species (EV-E/EV-F) and different subtypes (EV-E2/EV-E7, EV-E1/EV-E7, and EV-E3/EV-E6) in the same cattle herds, and revealed the enterovirus genetic exchange in intraspecies and interspecies recombination. These results provide an important update of enterovirus prevalence and epidemiological aspects and contribute to a better understanding of enterovirus genetic diversity, evolution, and pathogenesis.

Recombination is a pervasive process generating diversity in most viruses, especially in Enterovirus.Recombination between serotypes for EVs was first observed in poliovi rus in 1991 (10) and since then in a wide range of human enteroviruses.Among the members of Enterovirus A, the recombination occurred quite commonly in spite of lower recombination rates in EV-A71 in comparison with other subtypes in this species (11)(12)(13).Recombination rates in Enterovirus B are high, and a vast number of recombinant forms have been identified in this species (14)(15)(16).Recently, Luo et al. reported that BEV strains (GXNN2204 and GXGL2215) were, respectively, obtained from the genomic recombination of EV-E4 with EV-F3, and EV-E2 to EV-E4 (17).Recombination events in EVs are almost exclusively detected at the edges of the structural P1 region and within the nonstructural 5′UTR, P2, or P3 regions (14,18).Although recombination has been shown to occur across the whole genome, not all recombinants may be viable.Thus, there would be selection of recombinants with breakpoints at particular locations (19).Since recombination is usually demonstrated by showing phylogeny inconsistencies (20,21), the recombinant forms (RFs) can then be assigned based on clustering of the strains in the nonstructural region (16,22), and breakpoints can be detected by performing bootscanning (21,23).According to the criteria by the International Committee on Taxonomy of Viruses, members of species in the genus Enterovirus should have a high degree of amino acid identity (amino acid >70% in the polyprotein and amino acid >60% in P1) (1).In 1999, Oberste et al. found that the VP1 coding region carries major neutralization epitopes among the capsid proteins and is likely the optimal region for virus identification and molecular typing (24).Analysis of the VP1 sequences has been applied to genotype the most enteroviruses such as EV-A71 (11), CVA16 (25), CVA6 (26), and EV-E3 (27).Recently, Zell et al. suggested that enteroviruses were classified based on the amino acid sequence identities of the VP1 capsid protein, of which 70% to 85% were defined as heterologous types and >90% as homologous types (28).
In this study, we have isolated 27 strains of the BEV from cattle herds with diarrhea in various regions in China from 2012 to 2018, identified two novel EV-E subtypes in cattle herds in China, unveiled the coinfection of two enterovirus species and different subtypes in the same cattle herds, and revealed the enterovirus genetic exchange in intraspecies and interspecies recombination.These results provide an important update of the enterovirus prevalence and epidemiological aspects and contribute to a better understanding of enterovirus genetic diversity, evolution, and pathogenesis.

Isolation and characterization of the novel BEV strains
To isolate the potential BEV, the fecal samples collected from cattle herds with diarrhea in Jilin, Henan, and Ningxia provinces during 2012-2018 were processed.After infection, the Vero cells began to show cytopathic effects (CPEs) as early as 12 hours post-infection (hpi) and detached at 24-36 hpi compared to the uninfected cells (Fig. S1A and B).The isolates were further identified by RT-PCR with BEV-specific primers (5′UTR).After sequence analyses, the amplified fragments were turned out to contain the sequences of BEV.Information on these novel bovine virus strains is listed in Table 1.The TCID 50 titers for the isolated BEV strains were determined to range from 10 2 /0.1 to 10 5.5 /0.1 mL (Fig. S1C).

Decoding of the complete genome sequences for novel BEV strains
To decode the complete sequences for the novel BEV strains, a series of primers were designed individually for each strain and were used to amplify the genomic sequences (Table S1).After sequencing and joining the overlapped PCR-amplified fragments, the complete genomes for these viruses were obtained.Analyses showed that the lengths of complete genome sequences for these novel BEV strains varied from 7,433 to 7,470 nucleotides; the lengths of the 3′UTR were relatively diverse, ranging from 93 nucleotides to 136 nucleotides; and the lengths of the 5′UTR ranged from 810 nucleotides to 828 nucleotides with the relative conserved sequence.Additionally, the ORFs of these novel strains were highly conserved, ranging from 6,501 nucleotides to 6,531 nucleotides and encoding a polyprotein of 2,167-2,177 amino acids.Moreover, all strains shared a typical picornavirus genome organization.The newly identified viral strains and their GenBank accession numbers are listed in Table 2.

The new BEV strains shared high sequence identity with EV-E and EV-F reference strains
To determine the sequence identity of these novel BEV strains with the reference strains, a comprehensive analysis of the nucleotide and amino acid sequences of VP1, P1, and 2C + 3 CD was performed.As given in Table 3, five strains, namely, HeN-A2, NX-FY40, JL-HY61, JL-HY68, and JL-DH13, shared the nucleotide sequence identity of 80.3%-89.1% and amino acid sequence identity of 96.5%-97.5% with the VP1 of EV-E3 reference strain HY12 (KF748290), indicating that these five isolated strains are evolutionarily closer to EV-E3 strains.In Table 3, JL-JY12, JL-JY29, JL-JY41, and JL-GZL4 strains shared a relatively high nucleotide and amino acid sequence identity with the VP1 of the EV-E1 reference strain VG-5-27 (NC_001859), respectively, varying from 79.8% to 80.9% and 94.3% to 96.1%, indicating that these four isolated strains belong to the EV-E1 strains.HeN-B78, HeN-B84, and HeN-C4 strains from Henan province shared a relatively high nucleotide  S3).Taken together, the aforementioned results indicate these novel BEV strains are either close to EV-E or EV-F based on their sequence identities.

Phylogenetic analyses clustered the novel BEV strains to EV-E and EV-F species
To further confirm the phylogenetic status of these novel bovine enterovirus strains, phylogenetic analysis was performed based on the amino acid sequences of P1 and 2C + 3 CD using the maximum likelihood method.The representative sequences of all known types of EV-E and EV-F, and the representative sequences of all 12 enterovirus species and 3 rhinovirus within the Enterovirus genus were used as outgroup sequences.
The HY12 virus was used as a reference strain of EV-E3 for 2C + 3 CD phylogenetic tree analyses since it is the only EV-E3 virus strain with 2C + 3 CD sequences available.As shown in Fig. 1 and 2, 24 isolated strains were clustered in the same clade as EV-E, and 3 isolated strains, namely, HeN-A12, HeN-B62, and HeN-YR91, belong to the branch of EV-F strains.These results were congruent with the species and subtypes determined based on sequence identity, revealing the inconsistency of phylogenetic clustering between those of P1 and 2C + 3 CD.

Different subtypes and species of enterovirus detected in the same cattle herd
Further analysis on the enterovirus species showed that the HeN-B78, HeN-B84, and HeN-B62 strains from the same cattle herd were clustered to different enterovirus species, where HeN-B78 and HeN-B84 strains were grouped to EV-E and HeN-B62 clustered to EV-F.These results clearly demonstrated the mixed infections of EV-E and EV-F simultane ously in the same cattle herd.Similarly, HeN-A12 and HeN-A3 strains from another cattle The analysis models for ML methods for the P1 amino acid are LG + G + I. Viruses are marked with symbols as follows: • refers to the EV-E strains obtained in this study; ▲ refers to the EV-F strains obtained in this study; ◇ stands for the CEV-JL14 strain; ○ refers to the HY12 EV-E strain isolated from cattle; △ stands for the SD-S67 strain (an EV-F) isolated from goats.herd also belonged to different species, where HeN-A3 was grouped to EV-E and HeN-A12 was grouped to EV-F.These results demonstrated the unusual coinfection of EV-E and EV-F in the cattle herd.
It is interesting to note that different BEV subtype infections were also revealed in the same cattle herds.As shown in Fig. 1 and 3, JL-JY12, JL-JY29, and JL-JY7; HeN-A2 and HeN-A3; HeN-B78, and HeN-B84 and HeN-B81 from the same herds were clustered to different subtypes.These results demonstrated the complexity of BEV infection.

Recombination detected in the novel isolated bovine enterovirus strains
To characterize the genomic features and explore the origin of the isolated strains, recombination analyses were performed using SimPlot software and the RDP4 package • refers to the EV-E strains obtained in this study; ▲ refers to the EV-F strains obtained in this study; ◇ stands for the CEV-JL14 strain; ○ refers to the HY12 EV-E strain isolated from cattle; △ stands for the SD-S67 strain (an EV-F) isolated from goats.with multiple algorithms.The P values by the individual tools implemented in the RDP are shown in Table S4.The complete genome sequences of 24 isolated EV-E and 3 EV-F strains were selected as the query sequences for similarity and bootscan analysis with EV-E and EV-F reference sequences, respectively.As shown in Fig. 4A and B, the HeN-B62 strain resembled the EV-F7-AN12 strain in the sequence before the VP2 junction (nearby 1,160 bp), while it was more similar to the EV-F3-PS-87-Belfast strain than other strains for the fragment region (1,160-3,480 bp).Additionally, HeN-B62 contained an unidentified sequence that was apparently not related to the PS-87-Belfast strain in the fragment region (3,480-7,400 bp).Therefore, two recombination events possibly occurred in the HeN-B62 strain: one occurred between EV-F7-AN12 and EV-F3-PS-87-Belfast strains, with the breakpoint located in the VP2 gene, and the other occurred between EV-F3-PS-87-Belfast and unidentified sequences, with one breakpoint located in the 2A gene (Fig. 4B).The HeN-A12 (EV-F1) strain showed the highest degree of similarity to the EV-F1-SD-S67 strain in the P1 region.However, it contained an unidentified sequence that was apparently not related to the SD-S67 strain in the P2 and P3 regions (Fig. S2A).No apparent recombination was found in the HeN-YR91 strain with EV-F (data not shown).The HeN-A2 strain showed the highest degree of similarity to HY12 (EV-E3) in the P1 region.However, it contained an unidentified sequence that was apparently not related to the HY12 strain in the P2 and P3 regions (Fig. S2B).No apparent recombination events were found in other isolated EV-E strains (data not shown).

DISCUSSION
In this study, we have isolated 27 novel bovine enteroviruses from the major cattle-rais ing regions during 2012-2018 in China, explored the genetic diversity and recombina tion of the novel BEV strains, unveiled two new subtypes of EV-E designated as EV-E6 and EV-E7, and revealed the different BEV species or subtype infection on the same cattle herds.
BEV is the causative agent associated with an emerging disease in China, which leads to an important economic loss to the cattle industry.Based on virus genetic variability and molecular difference, BEV is divided into two species (genotype or serotype), EV-E and EV-F, where EV-E and EV-F are further divided into five subtypes (subgenotype) and seven subtypes, respectively (3).Although EV-E and EV-F infections were reported in China in the early 2010s (9,27), they remain largely unknown on many aspects such as the molecular difference, the virus genetic diversity and origin, the geographical distribution, and the transmission pattern.To unveil the underlying aspects of BEV infections, we performed BEV isolation using the representative specimens detected as BEV-positive by ELISA methods (29).Molecular characterization demonstrated that 27 BEV strains were successfully obtained.Sequencing and analyses of these BEV strains revealed that the majority of them are EV-E, suggesting that the predominant BEV infections are caused by EV-E.Interestingly, we found that different subtype strains within EV-E or EV-F were isolated from the same cattle herds; even EV-E (HeN-B78 and HeN-B84) and EV-F (HeN-B62) were both revealed in the same cattle herds.These results clearly demonstrated the coinfection of enterovirus species or subtypes.Since mutation and recombination extensively occurred among different BEV genotypes, the emergence and importation of new genotypes or subtypes of BEV should be continu ously monitored.The findings that the same subtypes were revealed in different regions or provinces indicate that BEV may be transmitted between different regions by live cattle distribution.More importantly, we found that some of these newly identified strains are likely the recombinants resulting from either enterovirus intraspecies or interspecies recombination, which enriches the understanding of underlying enterovirus pathogenicity and evolution.The amino acid sequence identity of VP1 is currently considered the standard for enterovirus classification (30).Recently, according to the suggestion by Zell et al., enteroviruses were classified based on the amino acid sequence identities of the VP1 capsid protein, of which 70% to 85% were defined as heterologous types and >90% as homologous types (28).Based on the criteria for enterovirus demarcation and our analyses on the available EV-E sequences, we proposed that the EV-E strains should be divided into seven subtypes.In addition to the five current reported subtypes EV-E1-EV-E5, obviously two novel subtypes for EV-E strains were discovered.We designated them as EV-E6 and EV-E7.This finding will broaden our current knowledge on enterovirus classification.
Genetic recombination plays important roles in the genetic variation and evolution of viruses.Although EV-A71 was mostly selected for recombination analysis in previous studies on Enterovirus, there were few studies on recombination analysis of bovine enterovirus (31)(32)(33).In the present study, the complete genome sequences of the novel BEV strains were used as the query sequences.Comprehensive recombination analyses revealed two recombination events that occurred possibly in the HeN-B62 strain.One recombination event occurred between EV-F7-AN12 and EV-F3-PS-87-Belfast strains, with the breakpoint located in the VP2 gene, and the other occurred between EV-F3-PS-87-Belfast and unidentified sequences, with the breakpoint located in the 2A gene, suggesting an intraspecies recombination.More importantly, interspecies recombination between EV-E and EV-F was found, where the EV-F7-AN12 strain is likely a recombinant from EV-E5 and EV-F1, and EV-E5-MexKSU/5 may be a recombinant from EV-F7 and EV-E1.It was also found that there was intraspecific recombination between the EV-F strains.We found that recombination breakpoints in the Enterovirus genome were frequently located in the nonstructural regions 5′UTR, P2, and P3 but rarely in the structural region P1, which is consistent with previous reports (34).The discovery that the nonstructural regions P2 and P3 act as the hot spots for recombination in BEV may be due to the higher homology of P2 and P3 regions than the structural region shared by BEV strains.Moreover, the inconsistent topology of phylogenetic trees suggested that it was not sufficient to genotype the recombinant viruses by VP1 alone.
The EV-F strains were first isolated in Northeast China in 2011.Since the discovery of the EV-E strain in 2012 in China, many BEV infections were reported.Our analyses on BEV strains demonstrated the EV-E strains as the predominant genotype in China.Since mutation and recombination extensively occurred among different BEV genotypes, the emergence and importation of new genotypes or subtypes of BEV should be continu ously monitored (9,27,30).

Sample collection and detection
A total of 2,721 fecal swabs collected from Jilin, Ningxia, and Henan were processed as previously described (27).All the samples were shipped on ice and then stored at −80°C.

TCID 50 titration
TCID 50 titration of the isolated strains was performed using 96-well plates following the procedure previously described (27).Briefly, the viruses were diluted at 10 × serial dilutions and used to infect the cells for each dilution.The cytopathic effects on the Vero cells were observed and counted at 48 hpi.TCID 50 was calculated following a standard procedure based on the Reed-Muench method (35).

RNA extraction and PCR
Viral RNA was extracted from the sample using the TRIzol reagent (Invitrogen) as previously described (27).cDNA was synthesized by reverse transcription using the Bio RT-cDNA kit (Invitrogen) following the manufacturer's instructions.The primers for amplification of the complete genome sequence were designed based on the variation of each strain (Table S1).cDNA ends were determined using the FirstChoice RLM-RACE kit (Thermos Fisher Scientific, Carlsbad, CA).PCR products were subjected to electro phoresis using 1.0% agarose gel and purified for DNA sequencing (Sangon Biotechnol ogy, Shanghai, China).The resulting sequences were analyzed using the Basic Local Alignment Search Tool (BLAST).The complete genome sequences were assembled using Lasergene software (version 7.0, DNASTAR, Madison, WI, USA).

Data collection and nucleotide sequence data sets
Genome sequences for enteroviruses from species E and F were retrieved from the GenBank database (to March 2023) (Table S2).Wrongly annotated sequences and coding sequences (CDSs) containing internal termination codons were discarded.A sequence alignment including the complete genomes of enterovirus strains was generated with Clustal W implemented in BioEdit, version 7.2.3,software (http://www.mbio.ncsu.edu/bioedit/bioedit.html)and split into three partitions corresponding to the VP1, P1, and 2C + 3 CD genomic regions.

Maximum likelihood phylogenetic analysis and recombination analysis
Maximum likelihood trees were generated with MEGA-7 software (36) using the best evolutionary model and 1,000 bootstrap replicates for the estimation of node consis tency (37).The merged sequence file (2C + 3 CD) was first aligned by MEGA-7 with the Clustal W method to generate nucleotide alignment.Phylogenetic trees were visualized by using the Interactive Tree of Life version 6 (https://itol.embl.de/).
Recombination analyses were performed using a recombination detection software package (RDP 4).A total of seven methods implemented in RDP4 were applied: RDP, Bootscan, GENECONV, 3Seq, Chimaera, MaxChi, and SiScan.Recombination detected by at least three of the seven methods with a P value cutoff 0.05 was considered true recombination.The similarity plot and bootscan analysis were performed by SimPlot 3.5.1 (JHK University, Baltimore, MD, USA).The pairwise similarity between all sequences in multiple sequence alignment was calculated with a 200-nt window moved along the sequence in 20-nt steps.

FIG 1
FIG 1 Phylogenetic analyses of P1 on the novel bovine enterovirus strains.The reference sequences include the representative sequences of all known EV-E and EV-F types and the representative sequences of all 12 enterovirus species and 3 rhinovirus species in the Enterovirus genus were used as outgroup sequences.The amino acid sequences of P1 were used to construct the phylogenetic tree using the maximum likelihood method with 1,000 bootstrap replications.Bootstrap values of >50 are shown at the nodes.The scale bar represents 10% nucleotide sequence divergence for maximum likelihood methods.

FIG 2
FIG 2 Phylogenetic analyses of 2C + 3 CD on novel bovine enterovirus strains.The reference sequences include the representative sequences of all known EV-E and EV-F types and the representative sequences of all 12 enterovirus species and 3 rhinovirus species in the Enterovirus genus were used as outgroup sequences.The amino acid sequences of 2C + 3 CD were used to construct the phylogenetic tree using the maximum likelihood method with 1000 bootstrap replications.Bootstrap values of >50 are shown at the nodes.The scale bar represents 5% nucleotide sequence divergence for maximum likelihood methods.The analysis models for ML methods for the 2C + 3 CD amino acid is LG + G + I. Viruses are marked with symbols as follows: • refers to the EV-E strains obtained in this study; ▲ refers to the EV-F

FIG 3
FIG 3 Phylogenetic analyses for BEV subtyping.Phylogenetic analysis was performed to subtype the BEV strain based on the VP1 amino acid sequences.The reference sequences including the representative sequences of all known EV-E and EV-F types and the representative sequences of all 12 enterovirus species and 3 rhinovirus species within the Enterovirus genus were used as outgroup sequences.A phylogenetic tree was generated using the maximum likelihood method with 1000 bootstrap replications.The scale bar represents 10% nucleotide sequence divergence for maximum likelihood methods.The analysis model for ML methods for the VP1 amino acid is LG + G. Viruses are marked with symbols as follows: red characters refer to the EV-E and EV-F strains obtained in this study, and white characters in the EV-E1 subtype refer to the EV-E strains obtained in this study; the green character stands for the strains isolated in the laboratory.The subtypes of the strains are shown in the figure.R-A stands for rhinovirus A, R-B stands for rhinovirus B, R-C stands for rhinovirus C.

FIG 4
FIG 4 Recombination analyses of the HeN-B62 strain with EV-F types.The genome sequence of HeN-B62 (A) was used as query sequences in the bootscan analysis.The structural map of the bovine enterovi rus genome is displayed below each panel.The possible crossover breakpoints were marked with a rectangular box with the enterovirus subtypes and represented by corresponding colors (B).Each point plotted is the percentage identity within a sliding window 200 nt wide centered on the position plotted, with a step size of 20 nt between points.

FIG 5
FIG 5 Recombination analyses between EV-E and EV-F types.The genome sequences of EV-F7-AN12 (A) and EV-E5-MexKSU/5 (B) were used as query sequences in the bootscan analysis.The structural map of the bovine enterovirus genome was displayed below each panel.The possible crossover breakpoints were marked with a rectangular box with the enterovirus subtypes and represented by corresponding colors.Each point plotted is the percentage identity within a sliding window 200 nt wide centered on the position plotted, with a step size of 20 nt between points.

TABLE 1
Molecular typing of enterovirus E and F isolates collected in China

Isolate Location EV type Years Isolation source Host Available sequence
EV-E strains, with the highest only 88.6% amino acid identity with EV-E2 strain PS 42 (DQ092792), indicating that those strains likely belong to a new subtype (Table

TABLE 2
The complete genome sequence of the newly identified virus strains

TABLE 3
Percent identity of nucleotide and amino acid sequences between isolated strains and enterovirus reference strains